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A Stirling engine configuration consisting of two cylinders, a regenerator and a sliding disc actuating 
mechanism (“swashplate”) is considered in this paper. A mathematical model, which combines funda¬ 
mental and empirical correlations, and principles of classical thermodynamics, mass and heat transfer 
accounting for variable heat transfer coefficients, is developed. The proposed model is then utilized to 
simulate numerically the system transient and steady state response under different operating and 
design conditions. A system global optimization for maximum performance in the search for optimal 
parameters that lead to maximum cycle efficiency is performed with low computational time. Appro¬ 
priate dimensionless groups are identified and the results presented in normalized charts for general 
application. The numerical results show that the two-way maximized system efficiency, ^ max ,max> occurs 
when two system characteristic parameters, the ratio between the total swept volume during the 
expansion, and the total swept volume, <p, and the ratio between the heat transfer area of the hot side 
heat exchanger and the total heat exchange area, y, are optimally selected, i.e., (<p,y) opt = (0.5,0.4). The 
two-way maximized cycle efficiency found with respect to the optimized parameters is sharp, in the 
sense that a 225% variation of the calculated efficiency values was observed within the range of tested 
configurations in this study, and “robust” (i.e., relatively insensitive) to the variation of several param¬ 
eters, thus stressing the importance to be considered in actual design. It is also found that the twice- 
maximized cycle efficiency and the total engine work output increase monotonically with the temper¬ 
ature of the hot source, Th. As a result, the model is expected to be a useful tool for simulation, design, 
and optimization of Stirling engines. 

© 2012 Elsevier Ltd. All rights reserved. 


1. Introduction 

Stirling engines are classified as thermal (or heat) engines 
within the theoretical thermodynamic framework. The Stirling 
engine thermodynamic cycle displays a theoretical thermal effi¬ 
ciency equal to the Carnot limit [1]. However, the heat transfer in 
this ideal limit must occur isothermally and reversibly, which 
demands infinite time, therefore, zero power is observed in the 
efficiency limit, l-T c /7h, in which Th is the temperature of the hot 
reservoir and T c the temperature of the cold reservoir. The 
mechanical difficulty to accomplish the ideal volume variation in 
the Stirling cycle is also responsible to separate the real from the 
ideal cycle. 
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A previous work by Curzon and Ahlborn [2] showed that the 
efficiency of a Carnot engine with a finite heat transfer time 
depends only upon the temperatures Th and T c , and is given by: 



As a result, the Stirling engine is a promising alternative, taking 
into account that it is an external combustion engine with possi¬ 
bilities for better combustion control and the potential to use 
multiple fuels. Therefore, with proper design, Stirling engines are 
expected to be less expensive and less polluting than Diesel engines 
and even gas turbines. 

The expectation of achieving thermal efficiencies close to the 
Carnot limit developed interest in the Stirling engine among 
researchers. Several authors have published theoretical and 
experimental studies on that possibility. Reader and Hooper [3] 
presented experimental measurements and analyses of Stirling 
engines. Uriel and Berchowitz [4] described analytically different 
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actuation mechanisms used in Stirling engines and presented 
mathematical models for the thermodynamic simulation of the 
cycle considering variable heat transfer coefficients and the pres¬ 
sure drop in tubular heat exchangers of circular cross-section. 

A finite time thermodynamic analysis of the Stirling cycle was 
presented by Ladas [5] with a sinusoidal volume variation and so 
called Beta configuration. The mass and energy balance equations 
were nondimensionalized and the maximum efficiency obtained 
was 14%. 

Finite time thermodynamics, endoreversible and simplified 
models have been used for Stirling engine optimization. For 
example: Blank and Wu [6] studied the power output and thermal 
efficiency of a finite time, optimized, extra-terrestrial, solar-radiant 
Stirling heat engine, in which the heat source and sink were assumed 
to have infinite heat-capacity rates, obtaining expressions for 
optimum power and efficiency at optimum power; heat pumps 
based on the reversed Stirling cycle were also studied [7] using 
a simplified model, which allowed a first optimization of real gas 
cycles, showing that efficiencies much higher than those achievable 
with an ideal gas, and similar to those of vapor-compression cycles 
can be obtained; a finite time thermodynamic optimization of 
a Stirling engine was performed by Popescu et al. [8], considering an 
endo- and exo-irreversible cycle, i.e., accounting for general irre¬ 
versibilities, finding optimum operating conditions leading to 
maximum power output in good agreement with experimental 
data; Erbay and Yavuz [9] studied the Stirling heat engine operating 
in a closed regenerative thermodynamic cycle finding the maximum 
power density and efficiency, in addition to the compression ratio at 
maximum power density; Bhattacharyya and Blank [10] reported 
major theoretical considerations concerning the design of an 
endoreversible Stirling cycle with ideal regeneration; Senft [ 11 ] used 
the classic Schmidt thermodynamic model (isothermal model) for 
Stirling engines and revisited the problem of identifying optimal 
engine geometry; Rogdakis et al. [12] performed the optimization of 
stable operation of the Stirling cycle in the free piston configuration, 
and for simplicity reasons used the Schmidt analysis; a thermoeco- 
nomic optimization of an irreversible Stirling cryogenic refrigerator 
cycle was presented by Tyagi et al. [13] finding that the effect of 
regenerative effectiveness is more than those of the other parame¬ 
ters on all the performance parameters of the cycle, and Mathieu 


et al. [14] performed the pre-sizing of the thermal storage and the 
solar receiver area of a thermodynamic micro power station and 
searched for possible optimum temperature levels. However, all 
such analyses were performed with steady state models. 

Two studies were found in the literature [15,16] that utilized 
dynamic models which included thermal losses to optimize the 
performance and design parameters of the Stirling engine, after 
experimental validation against data obtained from the General 
Motor GPU-3 Stirling engine prototype, showing significant engine 
efficiency increase with respect to the original configuration, i.e., 
from 39% to 51% increase. Although the study was dimensional and 
specific to the GPU-3 Stirling engine, it shows the importance of 
using dynamic models for engine optimization in order to produce 
more realistic results. For example, Karabulut [17] analyzed a free 
piston Stirling engine with a dynamic model, which made possible 
to show that the closed cycle performs a stable operation within 
a small range of the hot end temperature and damping coefficient 
of the piston motion that could be partially enlarged by inverting 
the engine into an open-cycle engine. 

The bibliographic review shows that there is a lack of studies on 
the thermodynamic optimization of Stirling engines with dynamic, 
dimensionless and more elaborated mathematical models. There¬ 
fore, this paper’s objective is to present a dimensionless dynamic 
mathematical model to simulate the thermodynamic behavior of 
a Stirling engine in the transient regime as a function of geometrical 
and operating parameters relevant to the engine design. Appro¬ 
priate dimensionless groups are introduced in order to present 
normalized simulation results for general engineering application. 
Acknowledging the finite availability of space in any engineering 
project, a constraint accounting for the total volume occupied by 
the engine is imposed. The mathematical model is then utilized to 
optimize some engine operating and design parameters for 
maximum cycle efficiency. 

2. Mathematical model 

Fig. 1 depicts a schematic diagram of the Stirling engine 
configuration considered in this work, consisting of: two cylinders, 
the hot and cold side heat exchangers, the regenerator and the 
sliding disc actuating mechanism (“swashplate”). One possible 


V d -2 



Fig. 1 . Schematic representation of the sliding disk mechanism (“swashplate”) and cylinders-regenerator unit consisting of a hot space, a cold space and a regenerator. 
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arrangement is a four-cylinder swashplate-drive double-acting 
engine (e.g., the Ford-Philips 4-215 engine [4]), which is the one 
considered in this work, i.e., each piston is responsible for the 
compression process in one side and the expansion process in the 
other side. 

The proposed thermodynamic model considers three control 
volumes: hot space (CV-1), cold space (CV-2), and the regenerator 
(CV-3). The model is dynamic, therefore in the analysis, tempera¬ 
tures in the three control volumes vary with respect to time, and 
volume in CV-1 and CV-2 also vary in time due to piston- 
displacement, but CV-3 undergoes constant volume processes 
since the regenerator volume is constant. The basic assumptions 
are: uniform temperatures in the control volumes, negligible 
pressure drops across the regenerators and heat exchangers, and 
the working fluid is treated as an ideal gas. The idea is to model the 
engine as an assembly of modules, like the one shown in Fig. 1, with 
identical thermodynamic performance. It is therefore sufficient to 
analyze one module and obtain the total engine power (or work) 
multiplying the power of one module by the engine’s number of 
modules. 

The pistons are driven by the swashplate mechanism illustrated 
in Fig. 1, resulting in a pure sinusoidal reciprocating motion having 
a phase difference between the adjacent pistons, given by the angle 
between pistons, a. Therefore, the total expansion volume, Vi, and 
the total compression volume, V 2 , vary in time, and are calculated 
by [4]: 


V 2 = V c + V A _ 2 + [1 + cos(wt)] (3) 

where Vh and V c represent the volume occupied by the hot and cold 
heat exchangers respectively, Vd-i and Vd -2 represent the dead 
volumes for expansion and compression, and V s _i and V s _2 the total 
swept volumes for expansion and compression for each cylinder. 

The phase angle between sets of cylinders, a , is shown in Fig. 2a, 
and the total number of cylinders, results directly from its value. For 
example, for the swashplate mechanism considered in this work, 
with a = 90°, the result is an engine with 4 cylinders and 4 modules 
as shown schematically in Fig. 2b. 

In Equations (2) and (3) the term (Vh + Vd_i) represents the total 
dead expansion volume and (V c + Vd- 2 ) the total dead compression 
volume. 

The time derivatives of the expansion and compression volumes 
for the sliding disk mechanism are given by: 

dVi Vc_i 

^ = — 2 1 oj sin(wt + a) (4) 


dV2 

dt 


V S -2 

2 


(o sin (cot) 



The volumes Vh and V c are functions of the heat exchanger 
geometry as follows: 



Vc 1 

^ + V d _! + 


[1 + cos (cot + a)] 



Vh = 


Ni • L\ -7T -D\ 





Fig. 2. (a) Illustration of the phase angle between piston units, a, and (b) The schematic representation of the components of one module in the double-acting swashplate Stirling 
engine considered in this work. 
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N 2 ‘ L 2 ' TU • 

4 ~ 



where N 1 , Li and Di are the number, length, and diameter of the 
tubes in the hot side heat exchanger, respectively, and Afe, I 2 and D 2 
are the number, length and diameter of the tubes in the cold side 
heat exchanger. 

Assuming uniform internal pressure, p, and neglecting the 
volume occupied by the gas inside the regenerator, the mass and 
energy conservation equations for the expansion and compression 
volumes state that: 


dm! 

dt 


m 



The temperature T R is found using the regenerator effectiveness, 
£, which is also a function of the flow direction and the tempera¬ 
tures 7i, I 2 and T m , in the following way: 

r R = a[T 2 + e(T m - T 2 )\ + b[T, +e(T m -h)} (17) 

The conditional temperatures, T c _i and T c _ 2 are, therefore given 
by: 

T c _i = a-T R + bT, (18) 

T c _2 = a-T 2 + f)-7p> (19) 

Equation (14) is rewritten, solving for dT m /dt: 


d (m-iT-i) 

L V-TT- 


Qi -P 


dV, 

dt 


+ c p mT c _i 



dr m 

dt 


Cp£ 

JTIrCr 


[ci(T2 — T m ) + b(T m 


h)] 


d 777 ] 

dt 


( 20 ) 


dm! dm 2 

-- H-- 

dt dt 


0 


The heat transfer rates in the expansion and compression 
(10) volumes are obtained as follows: 


Cv 


d(m 2 r 2 ) 

dt 




+ c p mT c _ 2 


( 11 ) 


where T c _! and T c _2 are conditional temperatures that depend on 
the gas flow direction inside the engine and on the heat exchange in 
the regenerator; rrq and T\ are the mass and temperature in the 
expansion volume, respectively; m 2 and T 2 are the equivalent 
quantities in the compression volume; c p and c v are the working 
fluid specific heat at constant pressure and at constant volume, 
respectively, and Q.! and Q .2 are the heat transfer rates in the hot 
and cold side heat exchangers, respectively. 

Combining the differential forms of the equations of state for each 
control volume, i.e., d(iri[T[)/dt = \/R(p(dV[/dt) + Vj(dp/dt)) for 
i = 1, 2, referring to CV-1 and CV-2, respectively, with Equations 
(9)—(11), it is possible to obtain differential equations to compute the 
internal pressure and the mass in the expansion volume as follows: 


Qi = 

^1^1 • (T h - Ti) 

(21) 

0.2 = 

^2^2 ' (Tc — T 2 ) 

(22) 


where h\ and /12 are the heat transfer convection coefficients, A! 
and A 2 the heat transfer areas in the expansion and compression 
volumes, respectively. The heat transfer areas A! and A 2 , are given 
by: 



= ■ l\ • TU • D! 

(23) 


= A^2 • L 2 * TU • D 2 

(24) 


The convection heat transfer coefficients vary during the engine 
operating cycle due to the variation in the gas mass flow rates and 
they are calculated using Colburn’s correlation [18] as follows: 


dp 

dt 


( 7 - 1 ) 


Qi 

T c -1 


+ 


0.2 
Tc-2 


+ py 


( 1 dVi 1 dV 2 \ 
\Jc- 1 dt + T c 2 dt ) 


Tc-, 


+ 


V 2 
Tc —2 


( 12 ) 


dm, = / V, \ dp / p \ dV, Q, 
dt \RyT c _J dt + \RT c _J dt c p r c _, 1 J 

where R is the gas constant and y the working fluid specific heats 
ratio c p /c v . 

The heat transfer rate between the regenerator matrix and the 
gas is computed as follows: 

Qr = m R c R ^ = c p [a(T 2 - T R ) + b(T R - T,)] ^ (14) 

where the coefficients a and b alternate between 0 and 1 according 
to the flow direction. T R is the temperature of the gas leaving the 
regenerator; T m the temperature of the regenerator matrix m R the 
matrix mass, and c R the regenerator specific heat. 

The values of the coefficients a and b are given by: 


h 


0.023 Rep /5 Pr 1/3 k 
D 


(25) 


where k is the gas thermal conductivity; Pr is the Prandtl number, 
and Reo is the Reynolds number based on the tube internal 
diameter. 

In order to proceed with the optimization of engine parameters 
for maximum efficiency, a volume constraint is introduced, to 
characterize the engine’s finite space availability: 

v* = Vc + V h + V d _! + V d _ 2 + V s _! + V s _ 2 (26) 

where V* is a reference volume, kept fixed during the optimization 
procedure. 

In order to nondimensionalize the mathematical model, the 
following scales are used: 


t* 

h* 


— d‘ 

0J 



m t Cv _ m t RT c 


V * 
D* 

T* 


Tc 



mcR 

m t c p 


(27) 


(a = 1) and (b = 0) for 
and 

(a = 0) and (b = 1) for < 0 


The dimensionless variables are therefore defined by 

(15) 

t = — t = — b = — 

L 1 1 t* r p* 

m i ,7 ^1 ,7 V 2 

(16) m i= - V, = v V 2 = v 


(28) 
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Additionally, the following dimensionless parameters for the 
engine are used: 


k 

Vi 


h D- 

h* 1 


Yi 

v * 


Li 


A A = A 

D* 1 A* 



(29) 


The index i, in Equations (28) and (29), refers to the points in 
s-1, s-2, d-1, d-2, c, h, 1, 2, m, R, c-1, c-2 in Fig. 1, representing 
the expansion swept volume, the compression swept volume, the 
dead expansion volume, the dead compression volume, cold side 
heat exchanger, hot side heat exchanger, expansion space, 
compression space, regenerator matrix, fluid in the regenerator, 
conditional temperature in the expansion space, and conditional 
temperature during compression, respectively, in the applicable 
variables. 

Substituting the dimensionless variables into Equations (12), 
(13) and (20) we obtain: 


4(V h + V c ) 

D i = --- L-, 

[y + (i -y)w]A 


(38) 


Using Eucken’s correlation, k = p(c p + 5.R/4) [19] for the gas 
thermal conductivity and Re d = mLD/(fiV) for the Reynolds 
number, as a function of the mass flow rate, tube length, tube 
diameter and volume occupied by the working fluid in the heat 
exchanger, it is possible to obtain expressions for h\ -A\ and 
h 2 -A 2 as functions of the dimensionless mass flow rate, drfq/dt, 
the gas properties and the dimensionless parameters previously 
defined: 


/i-[ "A-[ 


T V 3 (9 T _ 5) 2 / 3 

2739 



dm! 

dt 


4/5 


(39) 


dp 

dt 



+ 


h 2 • A 2 •^1 


T 
1 c 


-2 



f 1 dVi | 1 dU 2 \ 

\T C , dt t c 2 dt J 




(30) 


dm] 

dt 


vM + yp^--h v A v (t h -t, 


dt 


m 


dt M 


dt dt 


a(f 2 -Tm) + b(t m -fj 


1 


T T 


(31) 


C-1 


] dm! 
dt 


(32) 


The dimensionless temperatures f R , f c _! and f c _ 2 are defined 
as follows: 


t R = a 


To + e ( T m — t 


+ b 


T\ + si T m — T] 


T c -1 = o-T r + b-T^ 


(33) 


(34) 


b 2 -A 2 = 


T 1 / 3 ( 9t _ 5 ) 2/3 


27.39 


m t D2 J 


4(1 -yjA 


wDi 


1/5 


dmi 


dt 


4/5 


(40) 


The dimensionless network, VV, and heat transfer, Q. t , per cycle, 
are found by: 


d W ~ (dV\ dV 2 

= PI + 


dt 


dQi 

dt 


dt dt 


= h\ -A! (f h - ?! 


(41) 


(42) 


T c-2 = d m T 2 -\-b Tr (35) 

Several parameters, related to the volume distribution are 
needed to study and optimize the engine: a, <p , i//, x, and /3. This last 
parameter is introduced to study the influence in the engine 
performance of the relation between the reference volume and the 
piston diameter. Additionally, the following parameters related to 
the heat exchangers geometry are used: y, Z\ , Z 2 , and w. The 
nomenclature fully presents their definitions and mathematical 
expressions. 

It is possible to show that several parameters defined in the 
previous paragraph have the following relationships: 


w = 


y-(i -x) 

x(l -y) 


(36) 


(V h + V'c) = *(1 - <r) 


(37) 


Assigning a zero value to W and Q.! at the beginning of the cycle, 
the total dimensionless work W t and heat transfer Q.!_ t during the 
engine cycle are represented by the values of W and respec¬ 
tively, when t = 1 (end of the cycle). 

The cycle efficiency, rj , is therefore defined using W t ^nd Qi-t* 


Wt 

V = (y-\)‘ (43) 

0-1-t 

3. Numerical method and error analysis 

Equations (30)—(32) define a system of three ordinary differ¬ 
ential equations that was solved numerically with an adaptive time 
step fourth-fifth order Runge-Kutta method [20], together with 
Equations (33)—(35), (39), (40) in the dimensionless time interval 
(0 < t < 1). As a result, the thermodynamic properties of the gas 
during an operating cycle were obtained. Such dimensionless time 
interval is the time required to complete an engine cycle. 
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For each set of geometric and operating parameters, the system 
is solved iteratively for correct boundary conditions, i.e., requiring 
that the unknowns p, rrq and f m computed at the end of the cycle 
(t = 1) do not differ from their values at the beginning of the cycle 
since the working fluid is in the same thermodynamic state. This is 
done by computing a cycle relative error and satisfying the 
following convergence criterion [21]: 


^cycle 


U 


[ j+l 


U; 


U; 


< tol 


(44) 


where u = p, rrq or T m , tol is a convergence tolerance limit (e.g., 
0.01 or 1%), and j > 0 is the iteration counter. All three unknowns 
need to be tested according to the criterion of Equation (44) to 
ensure convergence. 

The process for finding correct boundary conditions starts from 
initial values estimated at the beginning of the cycle in the first 
iteration, which are replaced by the calculated values at the end of 
the first iteration, to start the second iteration and so on, until 
convergence is attained. In this work, the selected point to start the 
process computing the cycle relative error according to Equation 
(44) was the beginning of compression, with estimated initial 
values at j = 0, p 0 = 1, f m 0 = 1 and rrq 0 evaluated with the ideal 
gases equation of state using p 0 = 1, f m 0 = 1. 

In the simulation, H 2 was used as the working fluid, and its 
properties were taken from the technical literature [22]: 
Cp = 14.510 J/(kg K), R = 4.124 J/(kg K), n = 1,25 10" 5 kg/(m s). 
Additionally, the following values were used: T c = 300 K, 
w = 10 tt rad/s = 3000 rpm, a = 90° and D p = 0.1 m. 

Regarding the optimization of parameters, first a physical 
investigation is conducted to select the appropriate parameters to 
be optimized, which in the case of the present study was the pair 
(<P,y). Therefore two levels of optimization were carried out in this 
study for maximum system efficiency. The ranges of variation of 
each parameter to be optimized were: 0.25 < <p < 0.7 and 0.2 < y < 
0.7. The selected discretization in all ranges for the efficiency 
maximization was the coarsest set for which the optimal value of 
each parameter did not change as the sets became finer, while the 
relative error was kept below 1% in all cases [21]. 



0 


4. Results 

Initially, it was investigated the existence of optimal values for 
the volume parameters, a, cp and x, using the following values for 
the remaining parameters: p* = 10 bar, f3 = 2, \j/ = 0,9, A= 25 , 
y = 0.5, Z\=Z 2 = 400, f h = 3, M* = 8 and x = 0.5. Notice that once, 
(3 and D p , are given, the reference volume V* can be computed. The 
total mass of gas can also be obtained from the specification of p*, 
using Equation (25). 

In order to indicate the validity of the assumptions, numerical 
procedure and engine characteristics, the analysis starts by pre¬ 
senting a p x V t diagram, which is shown in Fig. 3a. The compres¬ 
sion and expansion processes are clearly devised in the graphic, 
showing that a positive network is achieved. For the engine 
configuration simulated in Fig. 3, the efficiency was rj = 24.25 %. 
Fig. 3b shows a p x 6 diagram, i.e., the behavior of the dimension¬ 
less indicated pressure with respect to crank angle during one 
engine cycle. Both graphs corroborate the expected trends for 
a Stirling engine, obtained with the simulations conducted with the 
mathematical model introduced in this paper. 

As the percentage of the reference volume occupied by the total 
swept volume increases, a (0,80, 0,85 and 0,90), it is observed an 
increase in the cycle efficiency as shown in Fig. 4. This analysis 
illustrates the variation of the cycle efficiency against the parameter 
cp , that represents the fraction of the expansion swept volume with 


Fig. 3. (a) The p x V t diagram of the double-acting swashplate Stirling engine, and (b) 
The p x 6 diagram of the double-acting swashplate Stirling engine. 
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Fig. 5. Behavior of the efficiency p with respect to <p, and varying x. 



Fig. 7. Maximized efficiency, 7? max , as a function of M*, and varying e. 


respect to the total swept volume. With an equal allocation of 
volume to the heat exchangers V h and V c , expressed by the 
parameter x = 0.5, it is observed a maximum efficiency corre¬ 
sponding to (p ~ 0.5, i.e., for equal swept volumes during 
compression and expansion. The optimal value of cp decreases 
slightly for smaller values of er. 

Setting cr = 0.9, the search for higher efficiencies continues by 
varying the parameter x, i.e., by varying the volume distribution 
between the heat exchangers. 

The simulation results are illustrated in Fig. 5. Higher effi¬ 
ciency values are obtained for larger values of x, i.e., when more 
volume is allocated to the hot end heat exchanger. It is important 
to note that x = ( V h /V h + \/ c ), therefore increasing x does not 
imply in an increase of the hot side dead volume, which is 
known to be a negatively affecting factor in the performance of 
Stirling engines. Increasing x means that more of the total 
available volume is allocated to the hot heat exchanger, therefore 
increasing the heat exchanger capability of collecting the avail¬ 
able heat supply. 

Once that the existence of a maximum efficiency corresponding 
to an optimal value of cp 0pt = 0.5 was demonstrated, the model was 
used to study the effect of the parameters Z\ and Z 2 , that represent 
the length to diameter aspect ratio of the heat exchanger tubes. The 
results, shown in Fig. 6, indicate that larger values of these 
parameters increase the engine efficiency. This trend is justified 


due to the increase in the flow velocity inside the heat exchanger 
tubes (which increases the heat transfer coefficient) and, the 
resulting diameter reduction, due to the fixed heat exchange area 
constraint, i.e., for y and A fixed. 

Fig. 7 illustrates the effect of the regenerator effectiveness, e, on 
the maximized efficiency which increases as e increases, which is 
expected since more energy is captured for work conversion. The 
graph also illustrates by means of the parameter M*, that the mass 
increase in the regenerator is unnecessary above M* = 1. 

Figs. (4)-(7) characterize the optimal value cp 0pt = 0.5 as 
“robust”, taking into consideration that it is insensitive to the 
variation of design parameters. This is valuable information from 
the point of view of engine design. 

Fig. 8 illustrates a second geometric optimization opportunity 
for the cycle. It shows that the cycle efficiency can be maximized 
one more time, now with respect to the parameter y, which 
represents the percentage of heat transfer area of the hot side heat 
exchanger with respect to the total available heat transfer area. The 
graph also illustrates the variation of the parameter p*, which is 
proportional to the total mass of working fluid used by the engine. 

It is found that as the total mass of working fluid is reduced, the 
twice-maximized efficiency increases. It is observed that the 
optimal value y op t = 0.4, is practically insensitive to variations ofp*. 
This finding is remarkable, as it characterizes the double optimi¬ 
zation as “robust”, i.e., independent of the main design parameters. 




Fig. 6. Behavior of the maximized efficiency, 7? max , with respect to Z 2 , and varying Z\. 


Fig. 8. Efficiency p as a function of y for variable p*. 
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Fig. 9. Efficiency p and work W t as functions of f h . 

Finally, the effect of varying the temperature of the hot source, 
7h, in the twice-maximized cycle efficiency is considered. The 
results in Fig. 9 indicate an increase in the efficiency with the 
increase in the temperature of the hot source. Fig. 9 also illustrates 
the engine total work, W t , that corresponds to the twice-maxi¬ 
mized efficiency, as a function of f h . 

5. Conclusions 

In this paper a mathematical model to simulate the operation of 
Stirling engines is introduced. A thermodynamic optimization of 
the cycle for maximum thermal efficiency is conducted. The 
numerical results show that there is an optimal engine geometry 
that maximizes the thermal efficiency. The results are reported 
using dimensionless variables, in normalized charts, which are 
general for the engine configuration considered, i.e., with a sliding 
disk mechanism (“swashplate”). 

The study identifies the location of the thermodynamic optima 
in terms of two important design parameters, cp and y, and their 
sensitivity to other engine parameters. The thermodynamic optima 
are sharp, in the sense that a 225% variation of the calculated effi¬ 
ciency values was observed within the range of tested configura¬ 
tions in this study, and “robust” (i.e., relatively insensitive) to the 
variation of several parameters, e.g., p* (or m t ), f h , M*, e , Z\, Z 2 , x and 
cr. It is therefore reasonable to state that this conclusion is of great 
importance for engine design. 
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Nomenclature 

A; total heat transfer area, Al + A2 

Ap. heat transfer area in the expansion volume 

A 2 : heat transfer area in the compression volume 

a, b: coefficients, Equations (15) and (16) 

c p : working fluid specific heat at constant pressure 

Cr: regenerator specific heat 

c v : working fluid specific heat at constant volume 

Dj: diameter of the tubes in the hot side heat exchanger 

D 2 : diameter of the tubes in the cold side heat exchanger 

hp heat transfer convection coefficient in the expansion volume 

hp heat transfer convection coefficient in the compression volume 

k: gas thermal conductivity 

Lp length of the tubes in the hot side heat exchanger 

Lp length of the tubes in the cold side heat exchanger 

mp mass in the expansion volume 

m 2 : mass in the compression volume 

m r. regenerator matrix mass 

mp total working fluid (gas) mass 

Np number of tubes in the hot side heat exchanger 

N 2 : number of tubes in the cold side heat exchanger 

Pr: Prandtl number 

p: internal pressure 

Q 4 : heat transfer rate in the hot side heat exchanger 
Qi 2 : heat transfer rate in the cold side heat exchanger 
Q r : heat transfer rate in the regenerator 
Q.; dimensionless heat transfer 
R: gas constant 

Rep: Reynolds number based on internal tube diameter 
t: time 

tol: cycle convergence tolerance 
Tp temperature in the expansion volume 
T 2 : temperature in the compression volume 
T c : temperature of the cold reservoir, 

T c _p conditional temperature that depends upon the gas flow direction and the heat 
exchange in the regenerator 

T c -2 : conditional temperature that depends upon the gas flow direction and the heat 
exchange in the regenerator 
Tp temperature of the hot reservoir 
T m : temperature of the regenerator matrix 
Tr: temperature of the gas leaving the regenerator 
u: general unknown 
Vp. total expansion volume 
V 2 : total compression volume 
V c : volume occupied by the cold heat exchanger 
Vd-P expansion dead volume 
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Vd -2 ■ compression dead volume 

VV volume occupied by the hot heat exchanger 

Vs_j; total swept volume for expansion 

V s _ 2 - total swept volume for compression 

w; diameter ratio D\ /D 2 

W: dimensionless cycle network 

x: ratio between the hot side heat exchanger volume, V h , and the total volume 
allocated to the heat exchangers, (V h + V c ) 
y: ratio between the heat transfer area of the hot side heat exchanger, A^ , and the 
total heat exchange area, (A^ + A 2 ) = A 
Zp. length to diameter ratio for the tubes of the hot side heat exchanger 
Z 2 : length to diameter ratio for the tubes of the cold side heat exchanger 

Greek letters 
a: phase angle 

(S; ratio between the reference volume V* and 7r.Dp 

y: specific heats ratio c p /c v 

e: regenerator effectiveness 

ecyde: cycle relative error 

p: cycle efficiency 

6: crank angle, cot x 180/7T 

a: fraction of the reference volume, V*, occupied by the total engine swept volume 
cp: ratio between the total swept volume during the expansion, V s _i, and the total 
swept volume, (V s _i + V^- 2 ) 


1 p: ratio between the volume allocated to the heat exchangers, (V h + V c ), and the 
total dead volume, (Vd-i + Vd -2 + + V c ) 

(a: angular speed (rad/s) 

Subscripts 

i used to represent the following: 
s—1: expansion swept volume 
s-2: compression swept volume 
d-1: dead expansion volume 
d-2: dead compression volume 
c; cold side heat exchanger 
h: hot side heat exchanger 
1: expansion space 
2 ; compression space 
m; regenerator matrix 
R: fluid in the regenerator 

c-1: conditional temperature in the expansion space 
c-2: conditional temperature in the compression space 
max: maximized (max, max = twice-maximized) 
t: total 

Superscripts 

*: reference value (scale) 
dimensionless variable 



